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1. INTRODUCTION 



This paper brings together the various scattered results on the 
Renewal Process model which are most relevant to applications. The 
early sections are expository in nature, describing the model and its 
applications, with particular emphasis on the Poisson process. Sub- 
sequent sections are devoted to general bounds and approximations for the 
renewal intensity, renewal function and forward recurrence distribution. 

A few special cases are worked out exactly and a numerical metiiod for 
the discrete case is given. Special approximations for the Wei bull 
case are included, and bounds on higher moments of the renewal counting 
process are given. Numerical illustrations are given whenever helpful. 

A theoretical development of the Renewal Process model is found in 
[Smith, 1958]. The entire subject matter, including some extentions, is 
fully discussed there and an extensive bibliography is provided. A good 
introduction to the subject is [Cox, 1962], where the fundamentals are 
treated in a very clear manner, for the case where inter-event times 
have a density. An excellent discussion of Renewal Theory is in [Feller, 1971], 
where a section is devoted to the two-sided case as well. We will not consider 
the more general theory based on two-sided distributions, as our applications 
call for the nonnegative case only. Some references to this generality are 
[Stone, 1965 and 1966]. Proofs of most results are not given here, but in- 
stead intuitive reasoning is provided to assist the reader's insight and 
comprehension. Finally, the problems of estimation of rates, tests for 
trends, and other issues of a statistical nature are left to [Cox and Lewis, 

1966] for resolution. This book gives an excellent discussion of the statistical 
problems which arise in practice. 



II. POINT PROCESS 

A point process is a series of epochs, t 3 • • • i n 

which describe the evolution of a process, usually stochastic in nature. 

This model is used to describe such processes as failure epochs of complex 
equipments, emissions of particles during decay or pulses along a nerve 
fibre for example. The renewal process is a special kind of point process. 
If the interevent times t^ » ^2 " t l ’ t 3 ‘ ^2’ "• ’ are stically 

independent random variables, and all except possibly the first have 
a common probability distribution, we have a renewal process. 

The simplifying assumptions of a renewal process make it a frequent 
choice as a model for certain point processes. One renewal process is most 
often assumed and its properties are particularly worth discussing, namely 
the Poisson process. It is the standard or benchmark to which other point 
processes are compared. 

III. THE POISSON PROCESS 

Consider a point process whose events occur at the epochs denoted by 
il'i-l^t^ ^.... If _ t-| , Xp - t^ - t-| , Xg = t.^ - t ? , . . . 
are statistically independent random variables each with the exponential 
di stri bution , 

Prob(X n s x} = 1 - e' Ax , x s 0 , 

then we call the process a (homogenous) Poisson process. The X are called 

n 

interevent times and the parameter A , the reciprocal of the mean (average) 
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time between events, is called the rate of the process, as it is the mean 
number of events occurring per unit time over any time interval. In fact, 
the actual number of events occurring per unit time, over an increasing time 
period, tends almost surely to A . 

There are many physical processes for which evidence suggests the 
Poisson process as an appropriate model. These include telephone calls 
made during a particular hour of the day, emissions from a radioactive 
source over a period during which the strength of the source remains substan- 
tially constant, the occurrence of nerve pulses along a nerve fibre, or 
the failures of complex equipment during its operating periods. See in 
particular [Drenick, 1960] for some general conditions under which the 
failures of complex equipment form approximately a Poisson process. Also, 
[Cox & Lewis, 1966, p. 251] show some data from several physical processes 
which exhibit the Poisson process assumptions. The following description 
is particularly suggestive of its wide applicability. 

Consider an arbitrary but fixed time horizon T , and a population 
of size n of individuals who may place a call during (0, T] . Suppose 
further that each individual decides to make a call with probability p , 
and that all individuals act independently, so that the number calling in 
(0, T] has the Binomial distribution with parameters n, p . Assume 
further that those calling select a time to call independently of others 
and uniformly from (0, T] . The number placing a call in, say (0, t] , 
with t <: T , has the Binomial distribution with parameters n , 
p(t/T) . Letting this random variable be N t , we note that 
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E(N t ) = At 

where 

A = np/T , 

and E denotes expectation. Now let n approach infinity and p 
approach zero such that A , the mean number of calls per unit time, 
remains constant. The limiting process is a Poisson process. The 
distribution of approaches the Poisson, mean At , and the times 
between successive calls are independent and have the exponential dis- 
tribution. One can imagine that a large population of potential callers, 
each one of whom has a very small probability of calling during (0, T] , 
and all of whom act independently, would approximate this limiting case 
quite well. The Poisson process is commonly referred to as the model 
for complete randomness in a point process. 

Another description of the Poisson process in terms of N^. will 
help to illustrate the model. Note that , the number of events 
during (0, t] , has a Poisson distribution, as this is the limit 
distribution for a Binomial with increasing n , decreasing p and con- 
stant mean np . Define N t ^ to be the number of events occurring 
in (t,t+h] , so that 



N t,t+h N t+h “ N t 



The following assumptions also characterize the Poisson process. 

(a) N^. has the Poisson distribution with parameter Ah , for 

all positive t and h . 
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is statistically independent of the number and position 



(b) N t, t+ h 

of the events in (0, t] , for all positive t and h . 

Some immediate consequences of (a) and (b) are: 

(1) Prob {N^ = n} = Prob (N^ = n } = e ^ ( A h ) n / n ! regardless 

of t , so that the process is stationary in time. 

2 

(2) Events occur singly since ProbiN^ ^ 2 } ~ (Ah) /2 near h = 0 . 

(3) In a short interval of length h , the probability of one event 
occurring is approximately Ah, and of no events, 1 - Ah . 

(4) The process' history prior to time t and its evolution after 

t are statistically independent. 

From (1) it follows that the first event occurs after h with 
A h 

probability e” , so that X-j , the time to the first event, has the 
exponential distribution. This reasoning applies to successive interevent 
times X by conditioning on the value of t , , so all interevent times 
have the exponential distribution. The independence of these times also 
follows from (4). 

To see that our definition is consistent with assumptions (a) and (b), 
we note that 

{N. < n} = {t > t) . 
t n 

This equivalence relates the counting variable N t to the epochs t and 
implicitly to the interevent times X since t = X, + .. + X . For 

the Poisson process, the X^ are iid with the exponential distribution, 
so t has tne Gamma distribution, yielding 
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Prob{N < n} = Prob{t > t} = )’ e Vt (At) k /k ! . 

r n k=0 

It follows that has the Poisson distribution, with mean At . The 

memoryless property of the exponential distribution is used to show (a) 
and (b) hold [Ross, 1972, p. 118]. 

The Poisson process enjoys many properties not shared by other point 
processes. These properties account in part for the wide use of the 
Poisson process in modelling problems. The counting variable N^. , 

the number of events in (t,t+h] , is of primary interest. It has the 
Poisson distribution with mean Ah and becomes nearly normal as Ah 
becomes large. By using a mean of Ah and standard deviation of /Ah , 
the normal approximation is quite reasonable when Ah ^ 20 . The 
counting variables corresponding to any set of disjoint intervals are 
statistically independent. This is particularly useful when accounting 
for the effects of events in different time periods, as the joint 
distribution of counts is easily specified. 

Any (non-random) number of independent Poisson processes can be super- 
imposed on a common time axis. The result is a Poisson process whose rate is 
the sum of the rates of the contributing processes. This follows directly 
from the reproductive property of the Poisson distribution under summation. 

It allows one to combine several processes without sacrificing any proper- 
ties. In fact, when many independent non-degenerate point processes, not 
necessarily Poisson, are combined, the pooled process usually tends to a 
Poisson process. This fact also accounts for the general acceptance of the 
Poisson process when modelling the output of several independent processes 
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of unknown origin. See [Cox & Smith, 1953, 1954] or [Khintchine, 1960] 
for a general discussion of the superposition of many processes and of 
its approach to the Poisson process. 

Another operation under which the Poisson process is invariant is 
the random deletion of events. Suppose that every event is subject to 
deletion with probability 1 - p , independent of other deletions. 

The remaining events form a Poisson process with rate pA , where X 
was the original rate. This result follows by either showing that the 
new interevent times are geometrical ly compounded exponentials and hence are 
still exponential, or by showing properties (a) and (b) hold by simple 
conditional probability arguments applied to the counting variable. This 
situation would arise if, for example, the original process could not be 
monitored directly, such as suicides from a bridge, or if one only 
recorded certain events, such as certain colored cars in a traffic 
stream; the recorded events would remain a Poisson process. As in the 
case of superposition, stationary point processes which are thinned 
this way become indistinguishable from Poisson Processes. By normalizing 
so the mean interevent time is one, the characteristic function for 
interevent times of the thinned process tends to 1 / ( 1 -0 ) , continous 
at the origin, as the probability p of event observation tends to 
zero. This corresponds to the exponential probability distribution with 
unit mean. Hence, series of events with substantial random loss should 
follow a Poisson process approximately. 

The epoch t at which the n^ 1 event occurs is also of interest, 
n 

As noted earlier, t has a gamma distribution with mean n/A and 
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2 

variance n/A . This is also equivalent to the fact that 2At n has 
the chi -squared distribution with 2n degrees of freedom, so standard 
tables can be used to find percentiles. This fact also provides a 
confidence limit or interval for A based on observation of t 

n 

When n becomes large, the central limit theorem shows that 
(At - n)//n has the standard normal distribution approximately. For 
improved accuracy, note that /At n is nearly normal with mean /n-1/4 
and variance 1/4. 

The conditional joint distribution of (t^,...,t ) , given the event 
{N t = n} , is the same as the distribution of the order statistics from 
n independent samples of a uniform distribution over [0,t] . This means 
that by conditioning on the number of events in an interval, the event 
positions in the interval have a relatively simple probabilistic law. 
Several statistical tests for deviations from the Poisson process are 
based on this [Epstein, I960]. Ross [1970, p. 18] uses this property 
to show that in a queue with an unlimited supply of servers and Poisson 
arrivals, the number of busy servers at time t is distributed Poisson 
wi th mean n(t) , 

t 

n(t) = A J o (1-G(x))dx , 

where A is the arrival rate and G is the distribution of service 
times. Since for large t the integral becomes E(S) , the mean 
service time, the steady state solution has the number of busy servers 
distributed Poisson with mean A E ( S ) . This model may arise in the 
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design of service systems where the number of servers is to be 
determined. The unlimited server case suggests the actual number 
needed for a low delay response. It may also arise where service is 
unlimited, such as for modelling the number recovering from industrial 
accidents. 

The local variables of a point process refer to those variables 
measuring the process at the time t . For the Poisson process, these 
variables have simple laws. The forward recurrence time, Y^ , is 
the time from t forward to the next event. Since Y t exceeds h 
when N^. = 0 , it follows that Y^ has the exponential distribution 

with parameter A . The backward recurrence time, , is the time 
since the last event prior to t , or t if none has occurred. By the 
same reasoning, has the exponential distribution with parameter 

A but is truncated at t . From property (b), Y t and are inde- 

pendent, so that the span, = Y^ + has, for t >> A - ^ , the 
gamma distribution with mean 2A~^ . This means that the interevent 
time which extends across the fixed time t is not a typical interevent 
time, as it has twice the mean of usual interevent times and a different 
distribution. This is a potential source of confusion and error in 
modelling problems, and it should be recognized at the outset. 

A Poisson process in multiple dimensions is called spatial in [Feller 
1968, p. 159]. Fie gives several examples of this case with data. 



- 9 - 



I V . RENEWAL PROCESSES 

When a Poisson process is generalized by allowing the interevent 
times to have a distribution distinct from the exponential, the result 
is a renewal process. To fix the notation, we use t , 0 <; t^ £ < . . . 

as the time of the n tn event, or renewal. The interevent times are 
X = t - t , ;> 0 , which are taken to be independent random variables. 

It is sometimes convenient to allow to have a distribution distinct 

from the others; however, for now we assume not, so that 

Prob{X n £ x} = F(x) , all n , 

and ProbtX > x) = F(x) . 

n v 

We also assume that F is a non-lattice distribution, except in the 
section on numerical methods for the discrete case. The second moment of 
F is assumed to be finite, with 

i 

A = J„ xdF(x) 

and = J o x^dF(x) - . 

In subsequent formulas, the product Aa , referred to as the 
coefficient of variation c , appears as a correction factor to the 
exponential case. It makes a convenient and reasonable measure of the 
variability of a renewal process. 

As before, we will use N^. as the number of events in (0,t] or 

= max{n| t p ^ t} . 
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For problems involving - N^. , we can refer to a renewal 

process whose first interevent time is Y^. , the forward recurrence 
time, as demonstrated later. 

From the central limit theorem applied to the sums t , and the 
equi val ence 

{N t < n) = {t n > t} , 

we have that, as t ■> 00 , 

N t - At d 

-> Normal (0,1 ) . 

Ao/At 

This means that, while the distribution of usually is intractable, 

a normal approximation for large t exists, and standard tables can be 
used to obtain percentiles. Note that the coefficient of variation appears 
as a correction factor to the standard deviation for the Poisson case. 

A quantity of considerable interest in modelling with a renewal 
process is the renewal function. 

M ( t ) = E(N t ) . 

This function, always non-negative and non-decreasing, appears in most 
formulas and is equivalent in information content to the distribution 
F . Both functions share common properties of differentiability, so that 
if F has a density, f , then M has a derivative, m , for which 

t 

M(t) = J 0 m(x)dx , 
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where m , the renewal intensity, may involve impulse terms if f 
does . 

A preliminary limiting result for the renewal function is that 

M(t) = At + o(t) 
as t 00 , 

where o(t) denotes a quantity tending to zero when divided by t . 

The only renewal process with a common interevent distribution for which 
rt(t) = At is the Poisson process. Under some regularity conditions 
(see [Smith, 1962]), the analogous result for the renewal intensity 
holds, namely 



m(t) = A + o(l) 
as t 00 , 

where o(l) denotes a quantity tending to zero. Both of these results 
are intuitively suggested by the definitions. 

For convenience, we will assume that the density f and renewal 
intensity m exists however the methods illustrated do not depend upon 
this. There are two lines of reasoning which can be identified that 
work well when modelling with a renewal process. The first involves 
conditional expectation methods and the second an interpretation of the 
renewal intensity. We will illustrate both methods, as they often 
complement each other. The "renewal equation" on which some treatments 
of renewal theory are based, is discussed first (see [Feller, 1941]). 
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To derive an equation for the renewal function, M , first condition 



on X.| , the time of the first renewal. We have that 



E(N t |X ] 




x £ t , 

X > t . 



Unconditioning by the distribution of X-j , we see that 

t 

M(t) = E(N. ) = F(t) + / i'l(t-x)f(x ) dx 

0 

and by differentiation 



t 

m(t) = f(t) + J m(t-x)f(x)dx . 



( 1 ) 



(la) 



We will refer to this as a renewal equation; several other equations 
which follow appear different but are equivalent. By taking the Laplace 
transform of this equation, (transforms being denoted by an asterisk*) 
we have 



or equivalently 



M* (s) 



f*(s) 

s(T-f*(s7T 



m* (s) 



f*(s) 

i~f*Ts7 



where for any g, g(t) = 0 for t < 0 , g*(s) = / e” St g(t)dt. 

0 

This equation serves to show that the renewal function and the 
interevent distribution are equivalent in information content, and hence 
suggests that equations involving a renewal process can be expressed with 
the renewal function. As we shall have occasion to use results for the 
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case where has a distribution different from the other interevent 
times, we leave as an exercise for the interested reader to show that 

t 

M^(t) = F-|(t) + J F-j (t-x)m(x)dx , (lb 

where quantities for the general case are subscripted by 1. 

The conditional expectation approach can be used to obtain the 
distribution of forward recurrence time. As earlier, we take Y to be 

the time elapsed from t until the next event, or t p - t for n = + 1 . 

Defining 

G(y,t) = Prob{Y t <; y) , 



we have by conditioning on the time until the first renewal 



Prob {Yj. £ y 



After unconditioning, we have 



f G(y,t-x) x <; t , 
1 =x}=\l t < x <: t + y , 



t 



t + y < x . 



G(y,t) = f o G(y,t-x)f(x)dx + F(t+y) - F(t) , 



which after some manipulation using the renewal equation (1) and 
transforms takes the form 

t+y 

G(y,t) = / F(t+y-x)m(x)dx . (2) 

t 

To illustrate the alternative method of derivation mentioned above, 
observe that m(x)6x is, to first order in 6x , M(x+6x) - M(x) , the 
mean of N x x+fix . As 6x -> 0 , this mean is Prob{N x x+6x = 1} + o(6x) . 
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In essence, the probability of multiple renewals in (x,x+fix] is of order 

(<Sx) , so at most one renewal need be accounted for and the mean 

number to occur is asymptotically the probability of exactly one renewal 
in (x,x+6x] , as fix -* 0 . Hence, m(x)6x can be viewed as the pro- 
bability of a renewal "at x". Equation (la) can now be reasoned by seeing 

that m(t-x)fit is the (conditional) probability of some renewal occurring 
at t , given the first renewal occurs at about x. From the law of 

total probability, and by passage to the limit in the choice of sub- 

t 

intervals for x , we have that / m(t-x)f(x)dx-6t is the probability 
of some renewal other than the first occurring at about t , to first order 
in fit. Now as f(t)5t is the probability of the first renewal occurring 
at about t , to first order in fit , we have 

t 

m(t)fit = (f(t) + { m(t-x)f (x )dx )fit + o(6t) ; 

upon dividing by fit and letting fit -*• 0 , equation (la) follows. 

Equation (2) is obtained in a similar manner, by observing that 
F(t+y-x)m(x)fix , to first order, is the probability that the last renewal 
prior to t+y occurs at about x , and that this event occurs at some 
x beyond t exactly when y . The case in which no density for 

the interevent distribution exists is handled by using the Lebesgue - 
Stieltjes integral and by replacing m(x)dx by dM(x) . 

Obtaining values for the distribution G(y,t) depends on being able 
to obtain the renewal intensity. In a later section, some approximations 
are discussed. Here, the steady state result for large t is given: 
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(3) 



y _ 

G(y.t) -* F (y) = A J F(x )dx 

0 

as t -* °° , 

where F g is called the equilibrium distribution corresponding to F . 

This result can be obtained from equation (2) by seeing that, for large 
t , m(x) = A in the range of integration (t,t+y] . Replacing 
m(x) by A in (2), we see that, as t , 

t+y. 

G(y,t) = j F ( t+y-x ) Adx = F (y ) . 
t 

A rigorous limiting arguement is the Key Renewal theorem which is proven 
in [Smith, 1958]. 

When the renewal process is first observed at a time t from its 
actual origin, the first interevent time is . For many practical 
situations, one must assume t is large enough to use steady state 
results for determining the distribution of Y^ . If we let X-| have 
the equilibrium distribution F g for F , while all the remaining 
interevent times have distribution F , then we have that 

t 

M e (t) = F e ( t ) + J 0 F e (t-x)m(x)dx , 

where M g (t) is the mean number of renewals in (0,t] for this equilibrium 
renewal process. Taking the Laplace transform of this equation gives 

Mg(s) = F g (s) + F e (s)m (s) 

and since * * 

F (s) 4(1 - f (s)) , 

e s 
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we have 



or 

M e (t) = At . 

Thus for the renewal process which is observed in steady state, the renewal 
function is At . 

The backwards recurrence time distribution can be obtained in the 
same manner as in the forward recurrence case. In particular, the same 
limiting distribution is obtained; 

Prob(Z, z) -* F (z) 
t e ' 

as t -> 00 , 

where is the time since the last renewal prior to t . Note that the 

exponential distribution is its own equilibrium distribution, and is unique 
in this respect. Finally, we note that our previous use of the renewal 
intensity immediately yields F(t-z)m(z) as the density of in (0,t) , 

with probability mass F(t) at z = t . Integrating this density gives 

t _ 

1 = F(t) + / F(t-z)m(z)dz (4) 

0 

which is equivalent to the renewal equation (1). 

The forward and backward recurrence times sum to the span, S t , 
which is the length of the interevent interval intersecting t . As noted 
in the Poisson process case, the span does not have the interevent distri- 
bution F . The sampling of interevent intervals is biased toward the longer 
intervals, as these are more likely to intersect the fixed time t . If 
h(s,t) is the density of , then for s < t , the probability h(s,t) • 6s 
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is the probability of a renewal occurring at x and the subsequent 
interevent interval being of length s , summed over x in [t-s,t] , 

or t 

h(s,t) = / f(s)m(x)dx = f (s)[M(t)-M(t-s )] . 
t-s 

As t tends to infinity, h(s,t) becomes h(s) , where 

h(s) = Asf(s) , 

the steady-state density of span. This result is obtained by observing 

that M(t) ~ At for large t . The result is actually an application of 

Blackwell's theorem, [Blackwell, 1948] . The limiting result can also be 

obtained by arguing that in a set of n random intervals, X-|,...,X , 

the combined expected length of all those of size x such that | x -s [ < 6_s , 

n "2 

is approximately n • s • f(s) • 6s which normalized by £ X. is the 

i=l 1 

probability of a span being within 6s of s . As n becomes large, we 
have the limiting result derived above for the asymptotic density of span 
by the strong law of large numbers. 

The remaining paragraphs of this section are devoted to the description 
of examples of applications involving the renewal process model. 

An early application of the model is to the theory of counters, as in 
[Feller, 1948]. A Geiger-Mu’l ler counter registers the emission of particles 
from a decaying substance, however due to the counter's resolving time, not 
all events are recorded. The evolution of recorded events can be modelled 
as a renewal process from which an estimate of the original process parameters 
can be made. 
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Models for a reliability problem frequently involve a renewal process. 
See [Gnedenko et. al., 1969] for a treatment of reliability theory using 
renewal theory. The interevent intervals correspond to machine operating 
periods and the renewal epochs to machine failures. The tail distribution 
for forward recurrence time, ProbiY^. > y) , is the interval reliability as 
defined in [Barlow and Proschan, 1965, p. 8], A model for the single 
item failure and repair process is the alternating renewal process. In 
this model, the process is in one of two states, 0 and 1 , and alternates 
between them. Residence times in state 0(1) are identically distributed 
samples from a distribution Fq(F-|) , and all residence times form an 
independent set. The epochs at which the process enters state 0 , or 
state 1 , form two imbedded renewal processes. The interevent distri- 
bution for either imbedded process is simply the convolution of Fq with 
F-| . To find the availability function, A(t) , which is the probability 
that the process is in, say state 1 at time t , assuming it began in 
state 1 , we have that 

t 

A(t) = F-j (t) + J Q [F 1 (t-x)]m(x)dx . 

The reasoning is that either state 1 was not left in [0,t] or, after 
reentering state 1 at time x , the process remained in state 1 beyond 
time t , for some time x in [0,t] . Some bounds in special cases 
are given in [Butterworth , 1963]; exact solutions are usually not available. 
Approximations can be found by using one of the approximations to m(x) 
reported below. The steady state availability follows immediately by 
replacing m(x) by A , its relevant value in the integrand when t is 
large, to obtain p, / (y ^ + Pq) in the limit. Here p.. is the mean 
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residence time in state i. The interval availability, A(t,t+h) defined 
as the probability of continued residence in state i over the interval 
[t,t+h] for an alternating renewal process, has an analogous derivation 
leading to 

t 

A(t,t+h) = F^(t+h) + | [F-| (t+h-x)]m(x)dx . (5) 

The steady state value in this case is just [u i / ( m -| + Pq)] [F^ e (h)] , 
which is suggested by the form of the expression. 

The study of replacement policies for failing equipment may involve 
the renewal model. Several bounds on the renewal functions are deter- 
mined in [Barlow and Proschan, 1964] (and by alternate methods in later 
sections here) by comparison of certain replacement policies. In 
particular, for a block replacement policy in which items are replaced 
every b units of time, and at failure, the expected number of failures 
between preventive replacements is M(b) . Accurate estimates of this may 
require knowledge of the renewal function for relatively short times 
b. The approximations given below are particularly directed toward the 
behavior of the renewal function near the origin. 

Bartholomew [1959, 1963^] applies the renewal model to problems in 
labor turnover, especially in circumstances where steady state-results 
are not of interest. For example, the renewal intensity is interpreted 
as the labor turnover rate for new organizations, and the interevent 
distribution as a means of measuring stability. An approximation to the 
renewal intensity for its transition from the origin to its steady state 
value is necessary for predicting the turnover effect. Actual steady state 
doesn't occur within the period of model applicability. 
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The queueing theory literature relies on the renewal process model 
to describe the input process to many service systems. Marshall and 
Wolff [1971] use bounds on the renewal function to get bounds on the 
difference between customer average and time average measurements for 
some queues. A related field. Inventory Theory, has problems involving 
a renewal process in the evaluation of the s - S policy (see [Arrow, 

Karlin & Scarf, 1958]). 

The following sections will give some bounds and approximations with 
practical significance, and briefly review some related literature. 

V. THE RENEWAL FUNCTION 

As we stated in section IV, a quantity which continually arises in 
renewal theory, and in the analysis of many models in which renewal 
processes are imbedded, is the renewal function, the expected number of 
renewals which occur in the interval (0,t], considered as a 
function of t . We defined this to be M(t) for an ordinary renewal 
process, and showed that 

t 

M(t) = F(t) + J M(t-x)f(x)dx , t ;> 0 . (1) 

o 

Another important representati on of M(t) is obtained by noting 
that the counting random variable (N t +1) is a stopping time in the sense 
of Wald, and by looking at the time to the first event after t , that is 
t + , after taking expectations 

M(t) = Xt + Xy(t) - 1 , tsO, (6) 
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where y(t) = E[Y t ] . 

a ) Approximations for t he_ Renewa 1 Function 
For all but a few simple processes it is difficult to solve 
equation (1) for M(t). It is equally difficult to determine M(t) in 
other ways, such as from (6), but since the renewal function appears in 
so many expressions for other quantities of interest, it is important that 
we can at least find approximations to M(t) . Early results emphasized 
asymptotic theory for large t . Clearly M(t) ^ 0 and it is non- 
decreasing. The simplest non-trivial result, often called the elementary 
renewal theorem (see, for example [Cox, 1962]) is that for any renewal 
process with finite mean interevent time 1/A , 

T M(t) 

1 im -► A , 

t -> oo z 

so that A can be interpreted as the "rate" of the process. 

Clearly M ( t ) increases "on the average" at rate A . The next 
result, which can be stated in a number of forms, shows that M(t) 
approaches a linear function with slope A (recall we are assuming that 
F is non-lattice. When F is lattice, consider M ( t ) to be defined 
only at multiples of the lattice; then M(t) is asymptotically linear). 
This well-known result is known as Blackwell's Theorem, and it states 
that for any h > 0 , with the restrictions just discussed, 

lim [M ( t+h ) - M ( t ) ] = Ah . 

t ->• oo 

p 

A similar result is found in [Smith, 1957] . Recall that a is the 

2 2 2 

variance of the inter-event times, and c = A o , the coefficient 
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of variation squared. Then if either F is non-lattice, or if t 
is restricted to points on the lattice, then 



M(t) = Xt + + o(l) . 



Let 




( 7 ) 



In this paper A. (t) will be used to denote the i th approximation 
to M(t) . Equation (7) can be used as an approximate formula for M(t) . 
It can be quite accurate for "large" t , where "large" depends on the 
distribution F as well as its mean. For many distributions, for t = 4 
or 5 multiple of tne mean, (7) is quite accurate. However if F is highly 
skewed, convergence to the linear function can be very slow. 

Equation (7) lias certain advantages and disadvantages. First, it 
is linear, easy to compute, and one needs knowledge of only the mean and 
variance of the interevent time. To offset these advantages, it can be 
quite poor for approximating M(t) for small t , and it is not known 
whetner trie true M(t) lies above as below A,(t) for any given t . 

Let us write A-j ( t ) as Xt + k , where k = . Now substitute 

this approximation to M(t) into the right hand side of the renewal 
equation (1) . After a little algebra the expression simplifies to 
A 2 (t) , where 



This gives us a second approximation for M(t) . It is easy to compute 
and has the correct asymptotic behavior; however, we have lost the linear 



A 2 (t) = Xt + (l+k)F(t) - F e (t) . 



( 8 ) 
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property. We are assured that A ? (t) is non-negative. This last property 

t 

is seen by noting that (1+k) ;> 0 , and F (t) = A J F(u)du ^ At . 

0 

From equation (1) it is easy to derive the following expression for 
M(t) eitner directly or by transform methods, 

M(t) = At + f o F e (t-u)m(u)du - F 0 (t) . (9) 



Equation (9) is useful for determining another approximation to M(t) . 

Let us replace m(u) on the right hand side of (9) by A , an exact 
expression for the case of a Poisson Process, but in general an approxi- 
mation. Then 

t _ 

A 3 (t) = At + A J o F e (u)du - F 0 (t) (10) 



is an approximation to M(t) . Like A^(t) it is no longer linear, and 
it requires knowledge of F and not just its moments. It is easy to show 
that 

c 2 -l 

lim (A 2 (t)-At) = j 1 , 

t -> 00 



so that tne approximation is asymptotically correct. 

The next approximation is derived from an approximation to the 
renewal intensity ni(t) derived in [Bartholomew, 1963]. Although it does 
not appear that the integral of the approximation was intended to be used 
as an approximation to M , it appears to give very good results for 
small t , although it may be poor for large t . In the notation of 
this paper, Bartholomew shows that 



m(t) ft; f(t) + 



\ F( t; 

A F e (t) 
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where f is the density of F . After a little algebra, the integral of 
his approximation (which we call A^(t)) can be written 



t 

A 4 (t) = At + F(t) - A J o 



F(u) 2 ' 

F e (u) 



du . 



do 



For many distributions it is computationally more difficult to use 
(11) than (10) or (8) because of the quotient in the integral. Note 
that the denominator F g is itself an integral of F . 

Let 



K(t) = J 
J 0 



F(u) 2 ] 
' F "e m 



du 



f F (u) - F(u ) 2 

" --F eW 



Let 

A(t) = F e (t) - F(t) . 



Then 

K(t) = J o F(u )du + J o A(u) • du . 

Now if F g (t) > F(t) for some t , then A(t)^|-^ < A(t) . 

Similarly if F (t) < F(t) the same inequality holds. Note also that 
since A(t) = F(t) - F g (t) , 

°° 1 c 2 

J t A(u)du - \- x - . 



Thus unless F (t) - F(t) VtsO, then 



K H < %- 
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From (11 ) we see that 



A 3 (t) = At + F(t) - AK(t) , 



and so if F (t) \ F(t) for some t , then 

c 2 1 

lim [A 3 (t)-At] > --■- - 1 . 

t °o 



Tne approximations given so far are not restricted to any particular 
family of distributions. Papers have appeared where F is restricted to 
a particular class. Leadbetter (1963] proves the following result. Suppose 
we can write F as an absolutely convergent power series 



Define 



F(t) = 



oo 



l 

n=l 



(-D n_1 

1'X 1 +mn] 



c t 
n 



mn 



m > 0 . 



= C 



1 



D 1 C 1 



Then 




n-1 

y D.C . 
j=l J n "J • 



H(t) = l 



H) 



n-1 



D_t 



mn 



fTT+mn] n 



m > 0 



( 12 ) 



where the series converges absolutely for all t s 0 . 

Applications of (12) to the case when F has a Wei bull distribution 
are given in [Leadbetter, 1963] and [Smith and Leadbetter, 1963] . In the 
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cases treated, only the first four terms of (12) were required. It should 



be mentioned that they treat higher moments of and in fact plot the 

variance of for various values of the parameter in the Weibull dis- 

tribution. We do not pursue higher moments in this paper. 

Further work on renewal processes for which F has a Weibull 
distribution is to be found in [Lomnicki , 1966], He discusses both the 
distribution of the number of renewals and the renewal function. We quote 
only his approximations to the renewal function. 

The basic idea in the work of Lomnicki is to write the renewal function 
as a power series whose terms are easy to calculate from available tabu- 
lated functions, and where convergence is rapid so that only a few terms 
need to be calculated. 

Assume that the underlying distribution F is given by 



The functions D are tabulated as the tail distribution of a Poisson 
random variable with mean t . Now define 



F(t) + 1 - e 



t ^ 0 , 3 > 0 . 



Define D.(t) = e _t l , k = 1,2,3, ... 

r=k r ' 



Y(r) „m, 




where 3 is the Weibull parameter. Next define 



bg(s) = y(s) ,s=0,l,2, ... 
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Next define 



Mp_> 

y(p) 




k = 0, 1, 2, ... 



s = k , k+1 , 



Finally let 



«k ( k ) a^(k) , 



s 



s-1 



a k (s) = l a r ( s ) - l a r (s) , s = k+1, k+2,... 



r=k r=k 



Then the renewal function 




s-1 3 k=l 



The interested reader should consult the Lomnicki paper for details such as 
convergence properties. The paper shows numerical calculations for a number 
of interesting functions in a Wei bull renewal process. 

b) B_ou_nds_for the Renewal Function 

Instead of approximations of the type presented in a) it is often 
desirable to nave upper and lower bounds on the renewal function, and a 
number of such bounds are presented here. It is usually easier to find 
lower bounds. Most upper bounds apply only with restricted classes of 
distributions, although some general upper bounds are given. 

The simplest non-trivial lower bound on M ( t ) is obtained from 
equation (6). By noting that y(t) is the expectation of a non-negative 
random variable we have our first lower bound L-j(t) , where 



L.j will be used for the ith lower bound, and 1). for the ith upper bound. 



L, (t) = At - 1 . 



( 13 ) 
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A lower bound was derived in [Barlow and Proschan, 1965, p. 68] 
in an interesting way, by comparing block and age replacement policies in 
a reliability theory context. The same bound, however, is easily obtain- 
able from equation (9) directly. Since F^ is a monotone non-increasing 
function, equation (9) gives 

M(t) ^ At + F e (t)M(t) - F e (t) . (14) 

Thus 

M(t) s L 2 (t) = F X -fvv - 1 . (14a) 

As was pointed out earlier, F g (t) ^ At , and so l_ 2 (t) iO . It is 
clear also that since F g (t) ^1 , that l_ 2 (t) ^ L-|(t) . 

Equation (1) can be used to find a simple upper bound on M(t) . Since 
F is monotone non-decreasing it is easy to see from (1) that 

M(t) s U,(t) = . (15) 

F(t) 

In general this bound is very poor and is of little practical use. The 
search for a good upper bound for a general process is more difficult. 

Stone [1972] in a short interesting paper derives the following very useful 
upper bound which has the nice properties of being linear in t and 
requiring only the first two moments of the underlying distribution; 

M(t) < U 2 (t) = At + 3(l+c 2 ) . (16) 

(Stone in this and a number of other papers deals with renewal processes 
in which the random variables can be negative. Since our motivation is 
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taken from reliability, inventory, and queuing theory where the random 
variables are invariably non-negative we do not pursue this two-sided 
general ization. ) 

If one is willing to restrict the class of distribution functions, 
other interesting bounds can be obtained. Many classes of distributions 
are described in [Marshall and Proschan, 1970], a paper which studies 
the relationship between various classes. 

Let us assume that the interevent distribution has the property 



The left hand side of this inequality can be thought to represent the 
mean residual lifetime of a component which has been in service a time 
t and which has a random lifetime distributed as F with unconditional 
mean 1/A . 

In reliability theory terminology we say that F belongs to the 
class of NBUE distributions (New-Better-than-Used-in-Expectation) , since 
clearly the inequality says that a used component can never have a mean 
residual life longer than the mean life of a new component. Any distri- 
bution for which the inequality in (17) is reversed is called UBNE with 
the obvious connotation. Equation (17) can also be written as 




(17) 



F g (t) £ F(t) , all t £ 0 . 



(17a) 



Equation (1) can be re-written in the form 



t 



F(t) = f F(t-u)m(u)du . 

o 



( 18 ) 
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Now if we assume that F is NBUE(UBNE) then using (17a) and (18) in 
(9) shows that 



M ( t ) £ At . 

(*) 



( 19 ) 



In this paper acronyms (such as UBNE) in parentheses should be read 
together with the inequalities in parentheses. Equations (13) and 
(19) together show that for F an NBUE distribution 



and we have very simple linear bounds. If F is UBNE (13) can be 
improved to M(t) £ At . 

A stronger assumption on F is to say that F has increasing 
failure rate (IFR). For details of IFR or other classes the reader 
should consult [Barlow and Proschan, 1965] or[Marshall and Proschan, 
1972]. Here we shall say that F is IFR if log F(t) is concave in 
t . For a renewal process with F restricted to this class, Barlow and 
Proschan [1965, p. 71], show that 



This result is obtained by comparing different replacement policies. No 
direct way of reproducing this result has been found. 

The linear bounds in (2) are very useful in many applications. Some 
are discussed later in this paper. If each bound is iterated in the basic 
renewal equation (1) , the reader will find that 



At - 1 <; M ( t ) £ At , 



( 20 ) 
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L 3 (t) = At - F g (t) <; M(t) s At - F g (t) + F(t) = U 3 (t) . 

It is interesting to note that the convex combination of these bounds 
which has the correct asymptotic behavior is the approximation A 3 (t) 
in equation (8) . 

Continued iteration of either bound yields a sequence of bounds 
which monotonical ly converges on M(t) from below or above. From a 
computational viewpoint this is an attractive property. The question 
naturally arises, can one start with better bounds than those in (20)? 

Of course the upper bound in (20) holds only for NBUE distributions, 
but Stone's upper bound could be used to start the iterative procedure. 

Marshall [1973] studied the problem of improving the starting 
bounds in this iterative procedure. It is easy to show that for any 
constant, say x , if the function At + x is used to start an iterative 
procedure in (1), the procedure will eventually converge on M(t) . 

Our problem is to find two numbers, b^ and b u such that 

At + b f s M(t) <; At + b^ , 

where At + b^ is as large as possible, and At + b u is as small as 
possible, and such that when either one is iterated in the renewal 

4 - 

equation (1), an improved bound is obtained for all t . These bounds 
are called the "best" linear bounds. 

To illustrate the method we use the lower bound. A simple iteration 
of At + by in (1) shows that the inequality 

t By "improved", we mean a bound which is no worse than the one in question. 
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At + b x ^ Xt + (l+b y )F(t) - F e (t) 



must hold if the improvement is taking place. 
Thus 



b 



£ 



FeU) 

F(t) 



l , 



and we make b. as large as we can. Let A be the set of t for which 
F(t) > 0 . Thus the reader should see that the largest b and smallest 
b u are given respectively by 



and 



F e (t) 

Inf 

teA F(t) 



Fg (t ) 

b = Sup - - 1 
u teA F(t) 



Since F e (0") = F(O') = 1 , -1 ^ b £ ^ 0 . Also if F is NBUE , then 
F g (t) <; F(t) , and so b u = 0 . In this case (19) is indeed the "best" 
linear upper bound. For some distributions b u is not finite, in which 
case M(t) cannot be bounded for all t by a linear function. 



c ) E xam ples and C alc ulations . 

For the Poisson process (i.e., F is exponential) the solution of 
(1) is M(t) = At . In this section (1) is solved for a number of renewal 
processes and compared to the approximations and bounds discussed in V 
a) and b). 

The first case considered is when F has a uniform distribution over 
(0, 2/A) . For this case (1) can be solved for successive intervals of 
length 2/A , and we find 
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M(t) = j; 



j-1 



n=0 




2 (j-1 )/X <; t <: 2 j / A 



j = 1, 2, 3, ... 



Note that the uniform distribution is in the NBUE class, and it is easy 
to see that the best linear bounds are given by (20) . Table 1 shows 
some of the approximations and bounds for the case X = 2 for t up to 
10 mean lifetimes. 

The second case shown is when F has a gamma distribution with mean 

2 

1 and variance -5 . As an example with C > 1 a third case is shown 

with F a hyperexponential distribution with mean 1 and variance 2.5 . 

Calculations for the gamma and exponential cases are shown in tables 2 

and 3 respectively. The gamma distribution is NBUE with k = - .25 , 

b = - .5 and b u = 0 . The hyperexponential distribution is UBNE with 

k = .75 , b = 0 and b = 1.5 . 
a u 
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TABLE 2. M(t) and approximations for the gamma distribution with mean 1 and variance 1/2. 



0.008 0.081 0.153 0.158 0.158 0.103 0.263 0.158 
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VI. THE RE NEWA L INTENSITY AND ITS CONVOLUTIONS. 

The previous section gives approximations and bounds on the renewal 
function. In this section, some approximations for the renewal intensity 
are shown. Also some approximations to convolutions of tail distribution 
functions with the renewal intensity are given and some applications for 
these are discussed. 

The renewal intensity function m(t) gives the rate of renewals 
occuring at the instant t from the start of a renewal process. If many 
(say N) identical renewal processes are operating simultaneously, then 
N * m(t) is the rate of renewals for the composite process. This fact 
is at the basis of an application to labour turnover. A new organization 
of fixed size experiences a turnover in its labour force which has been 
observed to fit this renewal model. See [Bartholomew, 1959] for a complete 
discussion with some data. 

As mentioned in an earlier section, m(t) approaches under suitable 
regularity conditions, the rate A for the process, as t tends to 
infinity. Our approximation is meant to describe the transition from 
M(0) = f(0) to the final value of A . Perhaps the best approximation 
we know of is the one given in the reference [Bartholomew, 1963 a ] , 
which we call h^(t) , where 



h-| (t) 



= f(t) + 



F* F(u )du 
0 



( 21 ) 



As seen later in numerical illustrations, h^ is an excellent 
approximation to m , particularly for highly skew interevent distributions. 



38 - 



f 



which is exactly when approximations are most needed. This is explained 
in part by the observations given by Bartholomew which are repeated here. 
The renewal equation (4) derived in section IV is essentially 

t 

1 = F(t) + | m(t-z)F(z)dz , 



or equalently 



t 

1 = F(t)/f m(t-z)F(z)dz . 
o 



From the renewal equation (la) we also have 



t 

m(t) = f(t) + f m(t-u)f(u)du . 
o 



Upon combining these, this renewal equation may be written 



t t 

m(t) = f(t) + F(t) / m(t-u)f (u)du/f m(t-u)F(u)du . 

0 0 



(la) 



The approximation h-| then results by replacing the ratio 

t t 

/ m(t-u)f (u)du/{ m(t-u)F(u)du 
0 0 

by 

t t _ 

f f(u)du/f F(u)du . 

0 0 

Since m(t) = A a constant for the exponential case, it follows that 
h.| = m for a Poisson process. Further, in cases where F is highly 
skewed so that m(t-u) in more nearly constant in u where f(u) and 
F(u) have appreciable positive value, the factoring and cancellation 
of m(t-u) is not so unreasonable and the approximation is close. 
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Other reasons for the success of h^ are: 

(1) h ] (0) = m(0) 

(2) Lim h-|(t) = Lim m(t) = X 

t °o t o° 

(3) hj(0) = ra'(O) and h^'(0) = m"(0) 

where prime denotes differentiation with respect to the single variable t. 

An advantage of h-j is that knowledge of the interevent distribution 

F is necessary only up to t; that is, h-j ( t ) is independent of the tail 

of F beyond t . This can be of considerable advantage when F is skewed 

and must be estimated from data. In particular, methods depending on 

transforms are inherently sensitive to the tail of F , a function which 

may not be well known in the upper regions, particularly if data is not 

available. The only possible drawback to h-, is the appearance of 
t _ 

/ q F(u)du which might have to be obtained by numerical integration. Even 
this however, should not be difficult. 

Another approximation to m(t) to consider is h 2 (t) where 

h 2 (t) = f(t) + AF(t) . (22) 

This approximation results by replacing m(t-u) by X in the right- 
hand side of equation (la) . Its relative lack of accuracy when com- 
pared to h-j is somewhat compensated for by the simpler functional 
form, however the choice of which one to use is best left to those making 
the application. While h 2 can itself be substituted in the rignthand 
side of (la) to give an improved approximation, the appearance of 
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convolutions quickly limits the practable usefulness of this approach, 
particularly in view of the availability of h^ . 

Tables 4 & 5 shown below give an indication of the numerical 
accuracy involved with these approximations. 



t 


m(t) 


h-j (t) 


h 2 (t) 


m(t) 


h-j ( t ) 


h 2 (t) 


0.0 


1 .60 


1.60 


1 .60 


4.51 


4.51 


4.51 


0.2 


1.51 


1.51 


1.36 


4.11 


4.11 


2.24 


0.4 


1.44 


1.44 


1.21 


3.75 


3.77 


1.40 


0.8 


1.32 


1.32 


1.04 


3.16 


3.25 


.99 


1.2 


1.23 


1.24 


.98 


2.69 


2.87 


.93 


1 .6 


1.17 


1.18 


.95 


2.32 


2.60 


.93 


2.0 


1 .12 


1.14 


.95 


2.04 


2.39 


.93 


2.5 


1 .08 


1.11 


.95 


1.76 


2.18 


.94 


3.0 


1 .05 


1 .08 


.96 


1.56 


2.02 


.94 


5.0 


1.01 


1.03 


.98 


1.17 


1 .62 


.95 


10. 


1 .00 


1 .00 


1.00 


1.01 


1.25 


.97 


X 


= 1-0, 


c = 1.58 




X 


li 

O 

s* 

n 


= 3.54 



TABLE 4. Calculations using a Hyper-Exponential interevent distribution 
with the given rate A and coefficient of variation c. 
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t 


m 


(t) 


h 1 (t) 


h 2 (t) 


0.0 


0 


• 


0. 


0. 


.01 




.039 


.039 


.039 


.02 




.077 


.077 


.078 


.04 




.148 


.148 


.151 


0.1 




.330 


.331 


.345 


0.2 




.551 


.556 


.598 


0.4 




.798 


.818 


.910 


1.0 




.982 


1.03 


1.14 


2.0 


1 


.00 


1.02 


1 .05 


4.0 


1 


.00 


1.00 


1.00 






A = 


1.0 c = 


0.77 



TABLE 5. Calculations using an Erlang interevent distribution with 
shape and scale parameters = 2.0 . 
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Another renewal process quantity to be approximated is a con- 



volution of the form 



t 

c(t) = / G(t-u)m(u)du 
o 



where G(x) = 1 - G(x) is any tail distribution function and m(u) is 
the renewal intensity. This convolution appears in application formulas 
such as for the availability in an alternating renewal process, the 
mean number present in a GI/G/ «> queue, and the distribution of the 
local variables forward and backward recurrence time. Since the forward 
recurrence time has the widest appeal to applications generally, we 
will illustrate both proposed approximations in the context of an 
approximation to the forward recurrence time distribution. 

The Key Renewal Theorem assures that, since the interevent dis- 
tribution F is assumed to be non-lattice, the limiting value for the 
convolution in question is 

t 

lim { G(t-u)m(u)du = yA , 

t - °° 0 

where 00 _ 

y = f G(x)dx (<°°) . 

J o 

(Note: y is the mean for distribution G .) 

To obtain our first approximation to the convolution c(t) , we 
can simply replace m(u) by h-j(u) or l^u) , and use numerical 
integration to obtain values. Since h^ is generally a closer approxima- 
tion, we take 

t _ 

c(t) ^ J G(t-u)h-| (u)du (23) 
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as an approximation to c(t) . As we seen in the numerical illustra- 
tions below, this can be a very close approximation when used for the 
forward recurrence time distribution. 

Another approximation to c(t) is easily obtained as follows. By 
considering G to be the distribution of ^ -j n the renewal process, 
and M-| to be the corresponding renewal function, we have the following 
version of equation (lb) 

t _ 

c(t) = / G(t-u)m(u)du = G(t) + M(t) - M,(t) . (24) 

0 l 

By replacing M(t) by A£(t) (see equation (8)) and M^(t) by an 
analogous approximation for the general case, we have 

M(t) ~ At + (l+k)F(t) - F e (t) 
and 

M-| (t) ~ At - F e (t) + G(t) + [l+k-Ap]F(t) , 
which upon substitution into equation (24) yield 

c(t) = ApF(t) (25) 



as an approximation to c(t) . 

Applying these approximations to the distribution of forward recur- 
rence time, , we have 

Prob(Y. > y) = F(t+y) + / F(t+y-u){f(u) + }du 

u 0 u __ 

f F(x)dx 
o 
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as the first approximation. Values are obtained by using numerical 
integration. 

To apply the second approximation to the forward recurrence time 
distribution, we can write 

t 

Prob{Y. > y} = F(t+y) + F(y) f G(t-u)m(u)du 
l o 

where 

G(x) = F(x+y)/F(y) , 

That is, G is the conditional distribution of residual life beyond y 
In this case, y is the mean residual life, or 

oo 

„ - /" ELXJ2). dx . . 

0 F(y) F(y) 

Now, substitution of equation (24) into the above equation yields 

Prob{Y t > y} ~ F ( t+y ) + F g (y)F(t) . 

The following table illustrates the two approximations in a 
particular case. 



(26b) 
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TABLE 6 





Hyper-Exponential 


2-Erlang 




Interevent Times 


Interevent limes 


y 


A = 1 .0 c = 1.58 


A = 1.0 c = .707 




t = 0.5 


t = 0.5 


0 . 


1 .0 


1.0 




1.0 


1.0 




1.0 


1 .0 


0.2 


0.765 


0.823 




0.765 


0.827 




0.804 


0.804 


0.4 


0.600 


0.653 




0.600 


0.657 




0.663 


0.629 


0.6 


0.483 


0.506 




0.483 


0.509 




0.558 


0.482 


0.8 


0.398 


0.385 




0.398 


0.388 




0.479 


0.363 


1 .0 


0.335 


0.289 


i 


0.335 


0.291 




0.419 


0.271 


1.5 


0.236 


0.135 




0.236 


0.135 




0.314 


0.125 


2.0 


0.179 


0.060 




0.179 


0.060 




0.247 


0.055 


2.9 


0.119 


0.013 




0.119 


0.013 




0.167 


0.012 



Forward Recurrence Time Distribution. Entries are Prob{Y. > y}; 
1st number is exact value, 2nd number using (26a), 3rd using (26b). 
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VII . CALCULATIONS IN DISCRJTE RENEWAL PROCESSES 

In many practical applications time is measured in discrete units, 
and the assumption of continuous distributions is an approximation. 
Renewal processes find application in inventory theory (for example 
see [Arrow, Karlin and Scarf, 1958] where the variable t is not 
"time" but discrete items of some commodity. In this section we present 
a simple computational method for solving renewal equations with discrete 
distributions and show some computational results. 

Let X. be the time between renewals (i-1) and i , and let 



F = Prob{X = n} , n = 0, 1, 2, ... , 
n n 



be the probability mass function, the same for all i . The ideas are 

first illustrated in obtaining the solution to the renewal function, 

equation (1). In the discrete case let be the expected number 

of renewals up to an including time n , and F = £ = Prob(X <; n} . 

n j=0 n 

Then equation (1) can be written 



M = F + y M .F. , n = 0, 1, 2, ... . 
n n n-j j 



( 1 ) 



For any fixed n , define the (n+1) x (n+1) lower triangular matrix 



0 0 

» 

fr 



0 



0 

'' f 



0 
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and let F and M be the (n+1) column vectors 



F = [F q , F 1 , F^, ... , F n ] , 
and 

M - [M 0 , M r M 2 , ... , M n ] 

respectively. Then equation (1) can be re-written as 

M = F + FM . 

Note that F is a non-negative matrix, every row sum is < 1 , and 
at least one row sums to strictly less than 1 (assuming E[X] > 0) . 
Therefore (I-F)”"* exists, where I is the identity matrix, and 

M = (I-F) _1 F . 

The inverse matrix of (I-F) has the structure 

G = 

and so is completely determined by a vector g = [g ...,g n ]. 

For small n G is probably quite easily calculated using a 
general matrix inversion routine. The primitive operator [T] in the 
APL language is extremely efficient. But for large n it is advan- 
tageous to make use of the special structure of F . 

Let us assume that (n+1) is a power of 2. If not we increase it 
to the next higher power. For example if (n+1) is 100 increase it 



g o, 0 ; 
9] 'V 

9 n • • • 



• 0 

s 0 
• 9 0 
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to 128 = 2 7 . Thus let us assume j is given by 2^ = (n+1), and 

write F. for the (n+1) x (n+1) matrix F , G. for the matrix G . 
J J 

If we partition F. into 



F. 

J 



F j-1 



H. , F. , 
J-l J-l 



, j = 1, 2, 



with Fq = Fq and Hq = Fi , then 



G 



j + 1 



G. 0 

J 

( G o H o G j) G j 



j - o, 1, 2, 



with G q = (l/l-f Q ) . 

Using an array manipulation language such as APL one can calculate 
large dimensioned matrices G in only a few iterations with no matrix 
inversion necessary. 

Figure 1 shows the renewal function for distribution bi no- 

mially with parameters 100 and .2 . Thus E(X.) = 20 , c = .2 
and k = -.48 . The approximation A^(t) = .05t-.48 is shown for 
comparison. A four-instruction program written in APL took only 1.8 
seconds on an IBM 370 model 158 to solve for M with n = 200 (10 mean 
times) . 

Having determined a method of solving (1) we have really solved 
a whole class of problems in renewal theory. In place of F let 
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6 



/ 





FIGURE 1 : Renewal .Function for a Binomi al Process wi th Parameters 

K)0 and 2. 
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us write a vector h - [h^, h-j , . . . , h ] . Now consider the renewal 
equation for the (n+1) - vector r , 

r = h + Fr , 

where F is the same matrix as before. Then 

r = Gh . 

By carefully choosing h a number of problems can be solved. Let 

h = [fg, f-j , ... , f ] . Then r is a vector of renewal intensities. 

Figure 2 shows the renewal intensity for a binomial renewal process 

with parameters 100 and .2 . As a further example consider the 

distribution of the excess random variable. Let Y be the time until 

n 

the first renewal after time n , and G(j,n) = Prob{Y n > j} , 

j = 0, 1, 2, ... . Define for a fixed j G. = [G(j, 0) , G(j, 1) , ... , 

G(j, n)] . If we let h = [F., Fj +1 , Fj + 2 , ... , F^ +n ] , then r will 

give us G. . In this way we can determine for a fixed j how the excess 
J 

distributions are converging with increasing n . 

Consider an alternating renewal process as described in Section IV. 

Let A n be the probability the system is operating at time n . Let the 

"up" times be distributed with distribution function U. = P [U . time ^ j] . 

0 u 

Let F. be the probability mass function of the sum of an up time and a down 
J 

time, and F the same as above. If h is taken to be the vector 
[U Q , U-| , ... , U ] then r will give the availabilities (assuming the 
system starts in the "up" state) at times 0, 1, ... , n . 
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FIGURE 2: Rene wal Intensity fo r a Binomial^ Process wjth Parameters 

100 and .2 . ' r 
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VIII. MISCELLANY. 

In this section we mention a few topics of special interest. 

These topics are higher moments of the renewal counting variable, some 
generalizations of the renewal model and superposition of several renewal 
processes . 

In some applications, higher moments of the renewal counting variable 
N.J. may be of interest. In particular, the variance of may be desired. 

Consider, for example, a situation in which events occur at times 
t i = i • a + e i , where e-j e^, ... , are I ID , normal (0, a) , as i 
runs over all intergers. This may be the arrival pattern of persons 
with appointments. This is clearly not a renewal process (unless a = 0) , 
however, when o is large, the point process appears, over fixed time 
periods, similar to a Poisson process. To distinguish this point process 
from the Poisson process or perhaps from other renewal processes, the 
variance of can be investigated. The variance of N t is also used 

in the analysis of the pooled output of several renewal processes, a 
topic mentioned below. 

Higher moments of are treated in [Cox, 1962, pp. 59-60] , 

t h 

where the Laplace transform methods are used. The r semi-invariant 
moment of is shown to be asymptotic to A^t + + 0(1) , (see 

also [Smith, 1959]) where A r is a function of the first r moments 
of the interevent distribution, and is a function of the first 

r + 1 moments. Iri particular, the variance of is (Xa) At +o(t) 

as t tends to 00 . 

2 

Since the quantity (Aa) appears as a factor times the variance 
of a Poisson Process with the same rate, it is reasonable in practical 
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as a measure 



investigations to use the coefficient of variation, Ao , 
of variation in the counting processes N t . Accordingly, when 
Ao < 1 (> 1) , we say the process is "under variance" ("over variance"), 
the boundary or bench mark case being the Poisson process for which 
Ao = 1 . 

Some bounds on higher moments are provided in [Barlow & Proschan, 
1964], where it is proven that 



E[U t (N t -l) ... (I4 t -k+l )] 

— kV" ' 




* U) 



(At) k 

kV 



when we have that F is NBUE (UBNE) respectively. It is also shown 
that the variance can be bounded by the renewal function, that is. 



V(N t ) s (>-) E(N t ) = M(t) 

when we have that F has increasing (decreasing) failure rate respectively. 

One generalization of the renewal process model discussed by Smith 
[1958] and Cox [1962] is the cumulative process. In this model, a random 
amount W. is contributed at each renewal epoch t. . We are primarily 
concerned with the net contribution made by time t . If represents 

this net contribution, we have 

N t 

Z t = I W i ’ 

1 i =1 1 

where the empty sum is given value zero. If the W. are identically 
equal to one, then = N. . If is the down time resulting from the 

i^ n replacement of a failed part, Z is the down time, cumulative to 
t . The W. may assume positive or negative values, if for example W. 
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is a monetary return associated with the i th event. The mean and 
variance of , when the W i are iid , and independent of N t , 
are given in [Cox, 1962, pp. 94] as 

E(Z t ) = M(t) * p w 
V(Z t ) = n(t) • ^ + V(N t )^ 

2 

wnere u , and o are the mean and variance respectively of the W. 

w w r J -\ 

distribution. Asymptotic normality of is also shown, with the 

asymptotic expansions for M ( t ) and V(N t ) replacing the exact 
values in the above expressions. 

In [Barlow & Proschan, 1964], the generalization in which the 
successive times between events may have different distributions, but 
must have a common mean, is discussed. Interevent times continue to be 
an independent set. The bound 

M(t) < (>) At 

is proven, when all interevent distributions have increasing (decreasing) 
failure rate respectively. These hypotheses can be relaxed to NBUE (UBNE) 
respectively. 

The superposition of several renewal processes provides a point 
process distinct from a renewal process, (except in the Poisson Process 
case). Such processes are observed when two or more arrival streams are 
pooled, for example, or when failures of a system are the failures of any 
component, and each component fails according to some renewal process. 
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The asymptotic result for many independent contributing processes is 
that, over intervals of fixed length, the composite process exhibits 
the properties of a Poisson process. The contributing processes may 
be non-identical, but the relative contribution to the total rate made 
by each contributing process must tend uniformily to zero. See 
[Drenick, 1960] for a complete statement of the necessary conditions. 

The superposition of only a few processes is treated in [Cox, 1962, 
pp. 71-77], wnere distributional results are emphasized, and in [Cox and 
Lewis, 1966, pp. 210-223] where a statistical analysis is given as well. 
See also [Cox and Smith, 1954] for a complete discussion. 
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